clear
set more off
capture postclose bskeep

use opioids_by_month

* construct variables used in analysis
keep if year>=2004
gen trend=12*(year-2004)+month
gen heroindr=(heroin_death/(population/100000))
gen opioidr=(opioid_deaths/(population/100000))
gen h_o_deathr=(h_o_deaths/(population/100000))
sort year
global mint=12*(2009-year[1])+1
global maxt=12*(2011-year[1])+12
local startnum = $mint
local endnum = $maxt
quietly summ trend if year==2004 & month==1
local firstnum = r(mean)
quietly summ trend if year==2014 & month==12
local lastnum = r(mean)

tempfile main bootsave 

qui save `main' , replace 

postfile bskeep trend sseh sseo sseb Fh Fo Fb using quadratic_spline_results_1, replace


forvalues b = $mint/$maxt { 
	display `b'
	use `main', replace 
	gen z=trend-`b'
	gen post=trend>=`b'
	gen zpost1=post*z
	gen zpost2=zpost1*zpost1
	gen zpre1=(1-post)*z
	gen zpre2=zpre1*zpre1
	reg heroindr zpre1 zpre2  zpost1 zpost2
	local sseh=e(rss)
	test (zpre1=zpost1)(zpre2=zpost2)
	local Fh = r(F)
	local ph = r(p)
	reg opioidr zpre1 zpre2  zpost1 zpost2
	local sseo=e(rss)
	test (zpre1=zpost1)(zpre2=zpost2)
	local Fo = r(F)
	local po = r(p)
	reg h_o_deathr zpre1 zpre2  zpost1 zpost2
	local sseb=e(rss)
	test (zpre1=zpost1)(zpre2=zpost2)
	local Fb = r(F)
	local pb = r(p)
	post bskeep (`b') (`sseh') (`sseo') (`sseb') (`Fh') (`Fo') (`Fb') 
}
postclose bskeep
clear
use quadratic_spline_results_1, clear
sort trend
list
desc
clear

use `main', replace
sort trend
merge 1:1 trend using quadratic_spline_results_1
egen minsseh=min(sseh)
egen minsseo=min(sseo)
egen minsseb=min(sseb)
gen mh=minsseh==sseh
gen mo=minsseo==sseo
gen mb=minsseb==sseb
egen th=max(mh*trend)
egen to=max(mo*trend)
egen tb=max(mb*trend)

gen z=trend-th
gen post=trend>=th
gen zpost1=post*z
gen zpost2=zpost1*zpost1
gen zpre1=(1-post)*z
gen zpre2=zpre1*zpre1
reg heroindr zpre1 zpre2 zpost1 zpost2
test (zpre1=zpost1) (zpre2=zpost2)
predict heriondrp

drop z post zpost* zpre*

gen z=trend-to
gen post=trend>=to
gen zpost1=post*z
gen zpost2=zpost1*zpost1
gen zpre1=(1-post)*z
gen zpre2=zpre1*zpre1
reg opioidr zpre1 zpre2 zpost1 zpost2
test (zpre1=zpost1) (zpre2=zpost2)
predict opioidrp

drop z post zpost* zpre*

gen z=trend-tb
gen post=trend>=tb
gen zpost1=post*z
gen zpost2=zpost1*zpost1
gen zpre1=(1-post)*z
gen zpre2=zpre1*zpre1
reg h_o_deathr zpre1 zpre2 zpost1 zpost2 
test (zpre1=zpost1) (zpre2=zpost2)
predict h_o_deathrp
drop z post zpost* zpre*

quietly summ trend if year==2010 & month==8
local dater = r(mean)
gen z=trend-`dater'
gen post=trend>=`dater'
gen zpost1=post*z
gen zpost2=zpost1*zpost1
gen zpre1=(1-post)*z
gen zpre2=zpre1*zpre1


** And figure out what Lambda is for significance table
gen startf = (`startnum' - `firstnum')/(`lastnum'-`firstnum' +1)
gen endf = (`endnum' - `firstnum')/(`lastnum'-`firstnum' +1)
gen lambda = endf*(1-startf)/( startf*(1-endf))
summ lambda

keep year month trend heroindr heriondrp sseh th lambda opioidr opioidrp sseo to h_o_deathr h_o_deathrp sseb tb Fh Fo Fb ph po pb
order year month trend heroindr heriondrp sseh th lambda opioidr opioidrp sseo to h_o_deathr h_o_deathrp sseb tb Fh Fo Fb ph po pb
outsheet using quadratic_spline_results_1.csv, comma replace


